options(stringsAsFactors = FALSE, scipen = 999)
try(Sys.setlocale("LC_CTYPE", "Chinese"), silent = TRUE)
script_arg <- grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE)
if (length(script_arg) == 1) setwd(dirname(normalizePath(sub("^--file=", "", script_arg))))

library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(patchwork)
library(broom)
library(purrr)
library(svglite)
library(ragg)
library(jsonlite)

dir.create("results", recursive = TRUE, showWarnings = FALSE)
dir.create("figures", recursive = TRUE, showWarnings = FALSE)

raw <- read_excel("Supplementary_File_1_Data.xlsx", sheet = "Data", skip = 1, .name_repair = "minimal")
stopifnot(ncol(raw) == 29, nrow(raw) == 173)

if (FALSE) {
names(raw) <- c(
  "study_id", "territory_text", "territory_code", "name", "age", "sex_code",
  "hypertension_code", "diabetes_code", "chd_code", "hyperlipidemia_code",
  "smoking_code", "alcohol_code", "education_le12_code", "nihss", "mrs",
  "visuospatial", "naming", "delayed_recall", "attention", "language",
  "abstraction", "orientation", "moca_raw", "legacy_psci", "legacy_severity",
  "laterality_raw", "site_raw", "corticality_raw", "multiplicity_raw"
)

aci_text <- "\u524D\u5FAA\u73AF"
pci_text <- "\u540E\u5FAA\u73AF"

clean <- raw %>%
  mutate(
    territory = case_when(
      territory_code == 1 | territory_text == aci_text ~ "Anterior circulation",
      territory_code == 2 | territory_text == pci_text ~ "Posterior circulation",
      TRUE ~ NA_character_
    ),
    territory = factor(territory, levels = c("Posterior circulation", "Anterior circulation")),
    territory_aci = if_else(territory == "Anterior circulation", 1, 0, missing = NA_real_),
    male = if_else(sex_code == 1, 1, 0, missing = NA_real_),
    hypertension = if_else(hypertension_code == 1, 1, 0, missing = NA_real_),
    diabetes = if_else(diabetes_code == 1, 1, 0, missing = NA_real_),
    chd = if_else(chd_code == 1, 1, 0, missing = NA_real_),
    hyperlipidemia = if_else(hyperlipidemia_code == 1, 1, 0, missing = NA_real_),
    smoking = if_else(smoking_code == 1, 1, 0, missing = NA_real_),
    heavy_alcohol = if_else(alcohol_code == 1, 1, 0, missing = NA_real_),
    education_le12 = if_else(education_le12_code == 1, 1, 0, missing = NA_real_),
    moca_adjusted = if_else(!is.na(moca_raw), pmin(30, moca_raw + education_le12), NA_real_),
    psci26 = if_else(!is.na(moca_adjusted), as.integer(moca_adjusted < 26), NA_integer_),
    psci22 = if_else(!is.na(moca_adjusted), as.integer(moca_adjusted < 22), NA_integer_),
    severity = case_when(
      is.na(moca_adjusted) ~ NA_character_,
      moca_adjusted >= 26 ~ "No PSCI",
      moca_adjusted >= 18 ~ "Mild",
      moca_adjusted >= 10 ~ "Moderate",
      TRUE ~ "Severe"
    ),
    severity = factor(severity, levels = c("No PSCI", "Mild", "Moderate", "Severe"), ordered = TRUE),
    severity_score = as.integer(severity) - 1L,
    corticality = case_when(
      corticality_raw %in% c("\u76AE\u5C42\u4E0B\u3001\u76AE\u5C42\u4E0B", "\u76AE\u5C42\u4E0B,\u76AE\u5C42\u4E0B") ~ "Subcortical",
      corticality_raw == "\u76AE\u5C42\u4E0B" ~ "Subcortical",
      corticality_raw == "\u76AE\u5C42" ~ "Cortical",
      grepl("\u76AE\u5C42", corticality_raw) & grepl("\u76AE\u5C42\u4E0B", corticality_raw) ~ "Cortical and subcortical",
      TRUE ~ as.character(corticality_raw)
    ),
    laterality = recode(as.character(laterality_raw),
                        "\u5DE6" = "Left", "\u53F3" = "Right", "\u53CC" = "Bilateral",
                        .default = as.character(laterality_raw)),
    multiplicity = recode(as.character(multiplicity_raw),
                           "\u5355\u53D1" = "Single", "\u591A\u53D1" = "Multiple",
                           .default = as.character(multiplicity_raw))
  )
}

# The distributed workbook is already deidentified and uses the standardized
# variable names documented on its Data Dictionary sheet. Recalculate every
# derived endpoint here so the published results remain independently checkable.
clean <- raw %>%
  mutate(
    territory = factor(territory, levels = c("Posterior circulation", "Anterior circulation")),
    territory_aci = if_else(territory == "Anterior circulation", 1, 0, missing = NA_real_),
    moca_recalculated = rowSums(across(c(visuospatial, naming, delayed_recall, attention, language, abstraction, orientation)), na.rm = FALSE),
    moca_adjusted_recalculated = if_else(!is.na(moca_raw), pmin(30, moca_raw + education_le12), NA_real_),
    psci26_recalculated = if_else(!is.na(moca_adjusted_recalculated), as.integer(moca_adjusted_recalculated < 26), NA_integer_),
    psci22_recalculated = if_else(!is.na(moca_adjusted_recalculated), as.integer(moca_adjusted_recalculated < 22), NA_integer_),
    severity_recalculated = case_when(
      is.na(moca_adjusted_recalculated) ~ NA_character_,
      moca_adjusted_recalculated >= 26 ~ "No PSCI",
      moca_adjusted_recalculated >= 18 ~ "Mild",
      moca_adjusted_recalculated >= 10 ~ "Moderate",
      TRUE ~ "Severe"
    )
  )

stopifnot(all(clean$moca_recalculated == clean$moca_raw, na.rm = TRUE))
stopifnot(all(clean$moca_adjusted_recalculated == clean$moca_adjusted, na.rm = TRUE))
stopifnot(all(clean$psci26_recalculated == clean$psci26, na.rm = TRUE))
stopifnot(all(clean$psci22_recalculated == clean$psci22, na.rm = TRUE))
stopifnot(all(clean$severity_recalculated == clean$severity, na.rm = TRUE))
clean <- clean %>% mutate(
  severity = factor(severity, levels = c("No PSCI", "Mild", "Moderate", "Severe"), ordered = TRUE),
  severity_score = as.integer(severity) - 1L
)

analytic <- clean %>% filter(!is.na(territory), !is.na(moca_adjusted))
lost <- clean %>% filter(is.na(moca_adjusted))

stopifnot(nrow(analytic) == 167, nrow(lost) == 6)
stopifnot(sum(analytic$territory == "Anterior circulation") == 91)
stopifnot(sum(analytic$territory == "Posterior circulation") == 76)

subtests <- c("visuospatial", "naming", "delayed_recall", "attention", "language", "abstraction", "orientation")
subtest_labels <- c(
  visuospatial = "Visuospatial/executive", naming = "Naming", delayed_recall = "Delayed recall",
  attention = "Attention", language = "Language", abstraction = "Abstraction", orientation = "Orientation"
)
subtest_max <- c(visuospatial = 5, naming = 3, delayed_recall = 5, attention = 6, language = 3, abstraction = 2, orientation = 6)

range_issues <- map_dfr(subtests, function(v) {
  bad <- which(!is.na(clean[[v]]) & (clean[[v]] < 0 | clean[[v]] > subtest_max[[v]]))
  tibble(variable = v, row = bad)
})
sum_mismatch <- analytic %>%
  mutate(component_sum = rowSums(across(all_of(subtests)))) %>%
  filter(component_sum != moca_raw) %>%
  select(participant_code, component_sum, moca_raw)

dup_nonmissing_ids <- clean %>% filter(!is.na(participant_code)) %>% count(participant_code) %>% filter(n > 1)

fmt_p <- function(p) {
  if (is.na(p)) return(NA_character_)
  if (p < 0.001) return("<0.001")
  sprintf("%.3f", p)
}
fmt_mean_sd <- function(x) sprintf("%.2f +/- %.2f", mean(x, na.rm = TRUE), sd(x, na.rm = TRUE))
fmt_med_iqr <- function(x, digits = 2) {
  q <- quantile(x, c(.25, .5, .75), na.rm = TRUE, names = FALSE, type = 2)
  paste0(format(round(q[2], digits), nsmall = ifelse(digits == 0, 0, digits)), " (",
         format(round(q[1], digits), nsmall = ifelse(digits == 0, 0, digits)), "-",
         format(round(q[3], digits), nsmall = ifelse(digits == 0, 0, digits)), ")")
}
fmt_n_pct <- function(x, denom) sprintf("%d (%.2f)", sum(x == 1, na.rm = TRUE), 100 * sum(x == 1, na.rm = TRUE) / denom)

group_values <- function(var) {
  list(
    anterior = analytic %>% filter(territory == "Anterior circulation") %>% pull(all_of(var)),
    posterior = analytic %>% filter(territory == "Posterior circulation") %>% pull(all_of(var))
  )
}

cat_p <- function(var) {
  tab <- table(analytic$territory, analytic[[var]])
  chi <- suppressWarnings(chisq.test(tab, correct = FALSE))
  if (any(chi$expected < 5)) fisher.test(tab)$p.value else chi$p.value
}

table1 <- bind_rows(
  tibble(Characteristic = "Age, years", Anterior = fmt_mean_sd(group_values("age")$anterior), Posterior = fmt_mean_sd(group_values("age")$posterior), `p value` = fmt_p(t.test(age ~ territory, data = analytic, var.equal = TRUE)$p.value), Test = "Independent-samples t-test"),
  tibble(Characteristic = "Admission NIHSS score", Anterior = fmt_med_iqr(group_values("nihss")$anterior, 0), Posterior = fmt_med_iqr(group_values("nihss")$posterior, 0), `p value` = fmt_p(wilcox.test(nihss ~ territory, data = analytic, exact = FALSE)$p.value), Test = "Mann-Whitney U test")
)
binary_t1 <- c(male = "Male sex", hypertension = "Hypertension", diabetes = "Diabetes mellitus", chd = "Coronary heart disease", hyperlipidemia = "Hyperlipidemia", smoking = "Smoking", heavy_alcohol = "Heavy alcohol use", education_le12 = "Education <=12 years")
for (v in names(binary_t1)) {
  av <- analytic %>% filter(territory == "Anterior circulation") %>% pull(all_of(v))
  pv <- analytic %>% filter(territory == "Posterior circulation") %>% pull(all_of(v))
  table1 <- bind_rows(table1, tibble(Characteristic = binary_t1[[v]], Anterior = fmt_n_pct(av, 91), Posterior = fmt_n_pct(pv, 76), `p value` = fmt_p(cat_p(v)), Test = ifelse(any(suppressWarnings(chisq.test(table(analytic$territory, analytic[[v]]), correct = FALSE))$expected < 5), "Fisher exact test", "Pearson chi-square test")))
}

or_from_glm <- function(model, term) {
  co <- coef(summary(model))[term, ]
  est <- unname(co["Estimate"])
  se <- unname(co["Std. Error"])
  tibble(term = term, OR = exp(est), CI_low = exp(est - 1.96 * se), CI_high = exp(est + 1.96 * se), p = unname(co["Pr(>|z|)"]))
}

primary_unadj <- glm(psci26 ~ territory_aci, data = analytic, family = binomial())
primary_adj <- glm(psci26 ~ territory_aci + age + male + education_le12 + nihss + hypertension + diabetes + hyperlipidemia, data = analytic, family = binomial())
sens_unadj <- glm(psci22 ~ territory_aci, data = analytic, family = binomial())
sens_adj <- glm(psci22 ~ territory_aci + age + male + education_le12 + nihss + hypertension + diabetes + hyperlipidemia, data = analytic, family = binomial())

primary_u <- or_from_glm(primary_unadj, "territory_aci")
primary_a <- or_from_glm(primary_adj, "territory_aci")
sens_u <- or_from_glm(sens_unadj, "territory_aci")
sens_a <- or_from_glm(sens_adj, "territory_aci")

severity_counts <- analytic %>% count(territory, severity, .drop = FALSE) %>% group_by(territory) %>% mutate(percent = 100 * n / sum(n)) %>% ungroup()
psci_counts <- analytic %>% count(territory, psci26) %>% group_by(territory) %>% mutate(percent = 100 * n / sum(n)) %>% ungroup()

tab2_rows <- list(
  tibble(Outcome = "PSCI, n (%)", Anterior = sprintf("%d (%.2f)", sum(analytic$psci26[analytic$territory == "Anterior circulation"]), 100 * mean(analytic$psci26[analytic$territory == "Anterior circulation"])), Posterior = sprintf("%d (%.2f)", sum(analytic$psci26[analytic$territory == "Posterior circulation"]), 100 * mean(analytic$psci26[analytic$territory == "Posterior circulation"])), `Effect / test` = sprintf("Unadjusted OR %.2f (95%% CI %.2f-%.2f)", primary_u$OR, primary_u$CI_low, primary_u$CI_high), `p value` = fmt_p(chisq.test(table(analytic$territory, analytic$psci26), correct = FALSE)$p.value)),
  tibble(Outcome = "No PSCI, n (%)", Anterior = sprintf("%d (%.2f)", sum(analytic$psci26[analytic$territory == "Anterior circulation"] == 0), 100 * mean(analytic$psci26[analytic$territory == "Anterior circulation"] == 0)), Posterior = sprintf("%d (%.2f)", sum(analytic$psci26[analytic$territory == "Posterior circulation"] == 0), 100 * mean(analytic$psci26[analytic$territory == "Posterior circulation"] == 0)), `Effect / test` = "", `p value` = "")
)
for (lev in c("Mild", "Moderate", "Severe")) {
  aa <- severity_counts %>% filter(territory == "Anterior circulation", severity == lev)
  pp <- severity_counts %>% filter(territory == "Posterior circulation", severity == lev)
  tab2_rows[[length(tab2_rows) + 1]] <- tibble(Outcome = paste0(lev, " PSCI, n (%)"), Anterior = sprintf("%d (%.2f)", aa$n, aa$percent), Posterior = sprintf("%d (%.2f)", pp$n, pp$percent), `Effect / test` = ifelse(lev == "Mild", "Ordinal comparison", ""), `p value` = ifelse(lev == "Mild", fmt_p(wilcox.test(severity_score ~ territory, data = analytic, exact = FALSE)$p.value), ""))
}
tab2_rows[[length(tab2_rows) + 1]] <- tibble(Outcome = "Education-adjusted MoCA score", Anterior = fmt_med_iqr(group_values("moca_adjusted")$anterior, 0), Posterior = fmt_med_iqr(group_values("moca_adjusted")$posterior, 0), `Effect / test` = "Mann-Whitney U test", `p value` = fmt_p(wilcox.test(moca_adjusted ~ territory, data = analytic, exact = FALSE)$p.value))
tab2_rows[[length(tab2_rows) + 1]] <- tibble(Outcome = "mRS score", Anterior = fmt_med_iqr(group_values("mrs")$anterior, 0), Posterior = fmt_med_iqr(group_values("mrs")$posterior, 0), `Effect / test` = "Mann-Whitney U test", `p value` = fmt_p(wilcox.test(mrs ~ territory, data = analytic, exact = FALSE)$p.value))
tab2_rows[[length(tab2_rows) + 1]] <- tibble(Outcome = "Anterior vs posterior circulation", Anterior = "Exposure", Posterior = "Reference", `Effect / test` = sprintf("Adjusted OR %.2f (95%% CI %.2f-%.2f)", primary_a$OR, primary_a$CI_low, primary_a$CI_high), `p value` = fmt_p(primary_a$p))
table2 <- bind_rows(tab2_rows)

score_results <- map_dfr(subtests, function(v) {
  a <- analytic %>% filter(territory == "Anterior circulation") %>% pull(all_of(v))
  p <- analytic %>% filter(territory == "Posterior circulation") %>% pull(all_of(v))
  tibble(variable = v, `MoCA subtest` = subtest_labels[[v]], `Anterior score, median (IQR)` = fmt_med_iqr(a, 0), `Posterior score, median (IQR)` = fmt_med_iqr(p, 0), score_p = wilcox.test(a, p, exact = FALSE)$p.value)
}) %>% mutate(score_q = p.adjust(score_p, method = "BH"))

nonfull_results <- map_dfr(subtests, function(v) {
  temp <- analytic %>% mutate(nonfull = as.integer(.data[[v]] < subtest_max[[v]]))
  tab <- table(temp$territory, temp$nonfull)
  chi <- suppressWarnings(chisq.test(tab, correct = FALSE))
  pval <- if (any(chi$expected < 5)) fisher.test(tab)$p.value else chi$p.value
  a <- temp %>% filter(territory == "Anterior circulation") %>% summarise(n = sum(nonfull), pct = 100 * mean(nonfull))
  b <- temp %>% filter(territory == "Posterior circulation") %>% summarise(n = sum(nonfull), pct = 100 * mean(nonfull))
  tibble(variable = v, `Anterior non-full, n (%)` = sprintf("%d (%.2f)", a$n, a$pct), `Posterior non-full, n (%)` = sprintf("%d (%.2f)", b$n, b$pct), nonfull_p = pval, anterior_pct = a$pct, posterior_pct = b$pct)
}) %>% mutate(nonfull_q = p.adjust(nonfull_p, method = "BH"))

table3 <- score_results %>% left_join(nonfull_results, by = "variable") %>%
  transmute(
    `MoCA subtest`, `Anterior score, median (IQR)`, `Posterior score, median (IQR)`,
    `Score raw p` = map_chr(score_p, fmt_p), `Score FDR q` = map_chr(score_q, fmt_p),
    `Anterior non-full, n (%)`, `Posterior non-full, n (%)`,
    `Non-full raw p` = map_chr(nonfull_p, fmt_p), `Non-full FDR q` = map_chr(nonfull_q, fmt_p)
  )

s1 <- tibble(
  Analysis = c("MoCA <22 prevalence", "Unadjusted association", "Adjusted association"),
  Anterior = c(sprintf("%d/91 (%.2f%%)", sum(analytic$psci22[analytic$territory == "Anterior circulation"]), 100 * mean(analytic$psci22[analytic$territory == "Anterior circulation"])), "Reference: exposure", "Reference: exposure"),
  Posterior = c(sprintf("%d/76 (%.2f%%)", sum(analytic$psci22[analytic$territory == "Posterior circulation"]), 100 * mean(analytic$psci22[analytic$territory == "Posterior circulation"])), "Reference group", "Reference group"),
  Estimate = c("Pearson chi-square", sprintf("OR %.2f (95%% CI %.2f-%.2f)", sens_u$OR, sens_u$CI_low, sens_u$CI_high), sprintf("aOR %.2f (95%% CI %.2f-%.2f)", sens_a$OR, sens_a$CI_low, sens_a$CI_high)),
  `p value` = c(fmt_p(chisq.test(table(analytic$territory, analytic$psci22), correct = FALSE)$p.value), fmt_p(sens_u$p), fmt_p(sens_a$p))
)

lesion_section <- function(var, label) {
  tmp <- analytic %>% mutate(category = as.character(.data[[var]])) %>% count(territory, category) %>% group_by(territory) %>% mutate(percent = 100 * n / sum(n)) %>% ungroup()
  ant <- tmp %>% filter(territory == "Anterior circulation") %>% select(category, n_ant = n, pct_ant = percent)
  post <- tmp %>% filter(territory == "Posterior circulation") %>% select(category, n_post = n, pct_post = percent)
  full_join(ant, post, by = "category") %>% mutate(across(c(n_ant, pct_ant, n_post, pct_post), ~replace_na(.x, 0))) %>%
    transmute(Characteristic = label, Category = category,
              `Anterior n (%)` = sprintf("%d (%.2f)", n_ant, pct_ant),
              `Posterior n (%)` = sprintf("%d (%.2f)", n_post, pct_post))
}
s2 <- bind_rows(
  lesion_section("laterality", "Laterality"), lesion_section("corticality", "Cortical involvement"),
  lesion_section("multiplicity", "Lesion multiplicity"), lesion_section("main_anatomical_site", "Main anatomical site")
)

s3 <- lost %>% transmute(
  participant_code, territory = as.character(territory), age,
  sex = if_else(male == 1, "Male", "Female"), nihss,
  education_le12 = if_else(education_le12 == 1, "Yes", "No"),
  hypertension = if_else(hypertension == 1, "Yes", "No"),
  diabetes = if_else(diabetes == 1, "Yes", "No"),
  coronary_heart_disease = if_else(chd == 1, "Yes", "No"),
  hyperlipidemia = if_else(hyperlipidemia == 1, "Yes", "No"),
  smoking = if_else(smoking == 1, "Yes", "No"), heavy_alcohol = if_else(heavy_alcohol == 1, "Yes", "No"),
  laterality, main_anatomical_site, corticality, multiplicity
)

model_table <- function(model) {
  tidy(model) %>% filter(term != "(Intercept)") %>%
    mutate(OR = exp(estimate), CI_low = exp(estimate - 1.96 * std.error), CI_high = exp(estimate + 1.96 * std.error),
           term = recode(term, territory_aci = "Anterior vs posterior circulation", age = "Age, per year", male = "Male sex", education_le12 = "Education <=12 years", nihss = "Admission NIHSS, per point", hypertension = "Hypertension", diabetes = "Diabetes mellitus", hyperlipidemia = "Hyperlipidemia")) %>%
    transmute(Variable = term, OR = round(OR, 3), `95% CI lower` = round(CI_low, 3), `95% CI upper` = round(CI_high, 3), `p value` = map_chr(p.value, fmt_p))
}
s4 <- model_table(primary_adj)

power_obs <- power.prop.test(n = mean(c(91, 76)), p1 = mean(analytic$psci26[analytic$territory == "Anterior circulation"]), p2 = mean(analytic$psci26[analytic$territory == "Posterior circulation"]), sig.level = 0.05, alternative = "two.sided")$power

deidentified <- clean %>%
  select(participant_code, territory, age, male, hypertension, diabetes, chd, hyperlipidemia, smoking, heavy_alcohol, education_le12, nihss, mrs,
         all_of(subtests), moca_raw, moca_adjusted, psci26, psci22, severity, laterality, main_anatomical_site, corticality, multiplicity)

data_dictionary <- tibble(
  Variable = names(deidentified),
  Description = c(
    "Deidentified participant code generated for this supplementary dataset", "Vascular territory", "Age in years", "Male sex (1=yes, 0=no)",
    "Hypertension (1=yes, 0=no)", "Diabetes mellitus (1=yes, 0=no)", "Coronary heart disease (1=yes, 0=no)", "Hyperlipidemia (1=yes, 0=no)",
    "Smoking history (1=yes, 0=no)", "Heavy alcohol use (1=yes, 0=no)", "Education <=12 years (1=yes, 0=no)", "Admission NIHSS score", "Three-month mRS score",
    "MoCA visuospatial/executive score (0-5)", "MoCA naming score (0-3)", "MoCA delayed recall score (0-5)", "MoCA attention score (0-6)",
    "MoCA language score (0-3)", "MoCA abstraction score (0-2)", "MoCA orientation score (0-6)", "Raw MoCA total", "Education-adjusted MoCA total, capped at 30",
    "Primary PSCI indicator: adjusted MoCA <26", "Sensitivity PSCI indicator: adjusted MoCA <22", "Programmatically derived PSCI severity", "Lesion laterality",
    "Main anatomical lesion site", "Cortical/subcortical classification", "Single/multiple acute lesions"
  )
)

write.csv(table1, "results/Table_1.csv", row.names = FALSE, fileEncoding = "UTF-8")
write.csv(table2, "results/Table_2.csv", row.names = FALSE, fileEncoding = "UTF-8")
write.csv(table3, "results/Table_3.csv", row.names = FALSE, fileEncoding = "UTF-8")
write.csv(s3, "results/Table_S1_Lost_to_Followup.csv", row.names = FALSE, fileEncoding = "UTF-8")
write.csv(s2, "results/Table_S2_Lesion_Distribution.csv", row.names = FALSE, fileEncoding = "UTF-8")
write.csv(s4, "results/Table_S3_Logistic_Regression.csv", row.names = FALSE, fileEncoding = "UTF-8")
write.csv(s1, "results/Table_S4_MoCA22_Sensitivity.csv", row.names = FALSE, fileEncoding = "UTF-8")
write.csv(deidentified, "results/Deidentified_Raw_Data.csv", row.names = FALSE, fileEncoding = "UTF-8", na = "")
write.csv(data_dictionary, "results/Data_Dictionary.csv", row.names = FALSE, fileEncoding = "UTF-8")

palette <- c("Anterior circulation" = "#3C5A97", "Posterior circulation" = "#D85845")
theme_pub <- function(base_size = 8.5) {
  theme_classic(base_size = base_size, base_family = "Arial") +
    theme(axis.line = element_line(linewidth = 0.4, colour = "#222222"), axis.ticks = element_line(linewidth = 0.4, colour = "#222222"),
          axis.text = element_text(colour = "#222222"), axis.title = element_text(colour = "#222222"),
          legend.title = element_blank(), legend.key.height = unit(3.2, "mm"), legend.key.width = unit(4.5, "mm"),
          legend.text = element_text(size = base_size - 0.6),
          plot.title = element_text(face = "bold", size = base_size + 1.2, colour = "#111111", margin = margin(b = 3)),
          plot.subtitle = element_text(size = base_size - 0.8, colour = "#444444", margin = margin(b = 4)),
          plot.tag = element_text(face = "bold", size = 12, colour = "#111111"),
          plot.tag.position = c(0.01, 0.99), panel.grid = element_blank(), plot.margin = margin(6, 7, 5, 7))
}
save_pub <- function(plot, stem, width_mm, height_mm) {
  w <- width_mm / 25.4; h <- height_mm / 25.4
  svglite(paste0(stem, ".svg"), width = w, height = h, system_fonts = list(sans = "Arial")); print(plot); dev.off()
  cairo_pdf(paste0(stem, ".pdf"), width = w, height = h, family = "Arial"); print(plot); dev.off()
  agg_tiff(paste0(stem, ".tiff"), width = w, height = h, units = "in", res = 600, background = "white"); print(plot); dev.off()
  agg_png(paste0(stem, ".png"), width = w, height = h, units = "in", res = 200, background = "white"); print(plot); dev.off()
}

# Figure contract: verified participant flow only; no unsupported screening or loss reasons.
flow_nodes <- tibble(
  xmin = c(0.25, 0.25, 0.06, 0.48, 0.425, 0.735), xmax = c(0.75, 0.75, 0.38, 0.94, 0.685, 0.995),
  ymin = c(0.82, 0.57, 0.29, 0.29, 0.05, 0.05), ymax = c(0.96, 0.71, 0.43, 0.43, 0.19, 0.19),
  label = c("173 patients enrolled\nFirst-ever MRI-confirmed acute ischemic stroke", "3-month follow-up assessment", "Incomplete follow-up\nn = 6; mRS and MoCA unavailable", "Included in analysis (n = 167)", "Anterior circulation\nn = 91", "Posterior circulation\nn = 76"),
  type = c("neutral", "neutral", "loss", "neutral", "aci", "pci")
)
flow_cols <- c(neutral = "#F5F6F7", loss = "#FCECEA", aci = "#E7EEF8", pci = "#FCEDEA")
fig1 <- ggplot() +
  geom_rect(data = flow_nodes, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = type, colour = type), linewidth = 0.65) +
  geom_text(data = flow_nodes, aes(x = (xmin+xmax)/2, y = (ymin+ymax)/2, label = label, colour = type), size = 3.2, lineheight = 1.08, fontface = "bold") +
  annotate("segment", x = .5, xend = .5, y = .82, yend = .71, arrow = arrow(length = unit(2, "mm")), linewidth = .45) +
  annotate("segment", x = .38, xend = .22, y = .57, yend = .43, arrow = arrow(length = unit(2, "mm")), linewidth = .45, colour = "#D85845") +
  annotate("segment", x = .62, xend = .71, y = .57, yend = .43, arrow = arrow(length = unit(2, "mm")), linewidth = .45) +
  annotate("segment", x = .635, xend = .555, y = .29, yend = .19, arrow = arrow(length = unit(2, "mm")), linewidth = .45, colour = "#315C94") +
  annotate("segment", x = .785, xend = .865, y = .29, yend = .19, arrow = arrow(length = unit(2, "mm")), linewidth = .45, colour = "#B94A3A") +
  scale_fill_manual(values = flow_cols) + scale_colour_manual(values = c(neutral="#30343B", loss="#B94A3A", aci="#315C94", pci="#B94A3A")) +
  coord_cartesian(xlim = c(0,1), ylim = c(0,1), clip = "off") + theme_void(base_family = "Arial") + theme(legend.position = "none", plot.margin = margin(7,7,7,7))
save_pub(fig1, "figures/Figure_1", 183, 105)

occ_plot <- psci_counts %>% mutate(Status = if_else(psci26 == 1, "PSCI", "No PSCI"), territory = factor(territory, levels = c("Posterior circulation","Anterior circulation"))) %>%
  ggplot(aes(x = territory, y = percent, fill = Status)) + geom_col(width = .62, colour = "white", linewidth = .45) +
  geom_text(aes(label = sprintf("%.1f%%", percent)), position = position_stack(vjust = .5), size = 2.9, colour = "#1F1F1F") + coord_flip() +
  scale_fill_manual(values = c("No PSCI"="#D6DCE5", "PSCI"="#D85845")) + labs(x = NULL, y = "Participants (%)", title = "PSCI occurrence", subtitle = "p = 0.012") + theme_pub() + theme(legend.position = "bottom", plot.subtitle = element_text(face = "italic"))
sev_plot <- severity_counts %>% mutate(territory = factor(territory, levels = c("Posterior circulation","Anterior circulation"))) %>%
  ggplot(aes(x = territory, y = percent, fill = severity)) + geom_col(width = .62, colour = "white", linewidth = .45) +
  geom_text(aes(label = ifelse(percent >= 7, sprintf("%.1f%%", percent), "")), position = position_stack(vjust = .5), size = 2.7, colour = "#1F1F1F") + coord_flip() +
  scale_fill_manual(values = c("No PSCI"="#D7DCE4", "Mild"="#9BB8CF", "Moderate"="#4F77A9", "Severe"="#D85845"), drop = FALSE) +
  labs(x = NULL, y = "Participants (%)", title = "PSCI severity distribution", subtitle = "ordinal p = 0.005") + theme_pub() + theme(legend.position = "bottom", plot.subtitle = element_text(face = "italic"))
fig2 <- (occ_plot | sev_plot) + plot_layout(guides = "keep") + plot_annotation(tag_levels = "A") &
  theme(legend.position = "bottom", plot.tag = element_text(face = "bold", size = 12), plot.tag.position = c(0.01, 0.99))
save_pub(fig2, "figures/Figure_2", 183, 94)

nf_long <- nonfull_results %>% select(variable, anterior_pct, posterior_pct, nonfull_q) %>%
  pivot_longer(c(anterior_pct, posterior_pct), names_to = "grp", values_to = "percent") %>%
  mutate(territory = recode(grp, anterior_pct = "Anterior circulation", posterior_pct = "Posterior circulation"),
         label = subtest_labels[variable], label = factor(label, levels = rev(unname(subtest_labels))))
nf_wide <- nonfull_results %>% mutate(label = factor(subtest_labels[variable], levels = rev(unname(subtest_labels))))
fig3 <- ggplot(nf_long, aes(x = percent, y = label, colour = territory, shape = territory)) +
  geom_segment(data = nf_wide, aes(x = posterior_pct, xend = anterior_pct, y = label, yend = label), inherit.aes = FALSE, linewidth = 1.25, colour = "#D0D4DA") +
  geom_point(size = 2.8, stroke = 0.7) +
  geom_text(data = nf_long %>% filter(territory == "Anterior circulation"), aes(label = sprintf("%.1f", percent)), size = 2.7, vjust = -0.85, show.legend = FALSE) +
  geom_text(data = nf_long %>% filter(territory == "Posterior circulation"), aes(label = sprintf("%.1f", percent)), size = 2.7, vjust = 1.55, show.legend = FALSE) +
  geom_text(data = nf_wide, aes(x = 104, y = label, label = paste0("FDR q ", map_chr(nonfull_q, fmt_p))), inherit.aes = FALSE, hjust = 0, size = 2.55, colour = "#555555") +
  scale_colour_manual(values = palette) + scale_shape_manual(values = c("Anterior circulation" = 16, "Posterior circulation" = 17)) +
  scale_x_continuous(limits = c(0,122), breaks = seq(0,100,20), expand = expansion(mult = c(.01,.01))) +
  labs(x = "Participants with a non-full subtest score (%)", y = NULL, title = "MoCA subtest profiles at 3 months") + theme_pub() +
  theme(legend.position = "top", legend.justification = "left", legend.box.spacing = unit(1.5, "mm"))
save_pub(fig3, "figures/Figure_3", 183, 108)

figs1 <- ggplot(analytic %>% mutate(territory = factor(territory, levels = c("Anterior circulation", "Posterior circulation"))), aes(x = territory, y = moca_adjusted, fill = territory)) +
  geom_boxplot(width = .48, outlier.shape = NA, alpha = .28, linewidth = .65) + geom_jitter(width = .11, size = 1.0, alpha = .55, colour = "#333333") +
  annotate("text", x = c(1,2), y = 31, label = c("n = 91", "n = 76"), size = 2.8, fontface = "bold") +
  scale_fill_manual(values = palette) + scale_y_continuous(limits = c(0,32), breaks = seq(0,30,5)) +
  labs(x = NULL, y = "Education-adjusted MoCA score", title = "Total MoCA scores at 3 months", subtitle = "p = 0.004") +
  theme_pub() + theme(legend.position = "none", axis.text.x = element_text(angle = 0), plot.subtitle = element_text(face = "italic"))
save_pub(figs1, "figures/Figure_S1", 89, 94)

forest <- tidy(primary_adj) %>% filter(term != "(Intercept)") %>%
  mutate(OR = exp(estimate), low = exp(estimate - 1.96*std.error), high = exp(estimate + 1.96*std.error),
         label = recode(term, territory_aci = "Anterior vs posterior", age = "Age, per year", male = "Male sex", education_le12 = "Education, 12 years or less", nihss = "Admission NIHSS", hypertension = "Hypertension", diabetes = "Diabetes", hyperlipidemia = "Hyperlipidemia"),
         label = factor(label, levels = rev(label)))
figs2 <- ggplot(forest, aes(x = OR, y = label)) + geom_vline(xintercept = 1, linetype = 2, colour = "#888888", linewidth = .4) +
  geom_errorbarh(aes(xmin = low, xmax = high), height = .14, linewidth = .65, colour = "#315C94") +
  geom_point(aes(colour = label == "Anterior vs posterior"), size = 2.4) +
  geom_text(aes(x = 18, label = sprintf("%.2f (%.2f-%.2f)", OR, low, high)), hjust = 0, size = 2.7, colour = "#333333") +
  scale_colour_manual(values = c(`TRUE` = "#D85845", `FALSE` = "#315C94"), guide = "none") +
  scale_x_log10(limits = c(0.2, 50), breaks = c(.25,.5,1,2,4,8,16), labels = c("0.25","0.5","1","2","4","8","16")) +
  labs(x = "Adjusted odds ratio (95% CI)", y = NULL, title = "Primary multivariable logistic regression") + theme_pub() +
  theme(axis.text.y = element_text(colour = "#222222"))
save_pub(figs2, "figures/Figure_S2", 150, 104)

summary_out <- list(
  data_checks = list(rows = nrow(clean), analytic = nrow(analytic), lost = nrow(lost), range_issues = nrow(range_issues), component_sum_mismatches = nrow(sum_mismatch), duplicated_nonmissing_ids = nrow(dup_nonmissing_ids)),
  primary = list(anterior_n = sum(analytic$psci26[analytic$territory == "Anterior circulation"]), posterior_n = sum(analytic$psci26[analytic$territory == "Posterior circulation"]), unadjusted_OR = primary_u, adjusted_OR = primary_a),
  sensitivity22 = list(anterior_n = sum(analytic$psci22[analytic$territory == "Anterior circulation"]), posterior_n = sum(analytic$psci22[analytic$territory == "Posterior circulation"]), unadjusted_OR = sens_u, adjusted_OR = sens_a),
  observed_effect_posthoc_power = power_obs,
  age = list(anterior = fmt_mean_sd(group_values("age")$anterior), posterior = fmt_mean_sd(group_values("age")$posterior), p = t.test(age ~ territory, data = analytic, var.equal = TRUE)$p.value),
  severity_p = wilcox.test(severity_score ~ territory, data = analytic, exact = FALSE)$p.value,
  moca_total_p = wilcox.test(moca_adjusted ~ territory, data = analytic, exact = FALSE)$p.value,
  mrs_p = wilcox.test(mrs ~ territory, data = analytic, exact = FALSE)$p.value,
  fdr_score = setNames(score_results$score_q, score_results$variable),
  fdr_nonfull = setNames(nonfull_results$nonfull_q, nonfull_results$variable)
)
write_json(summary_out, "results/analysis_summary.json", pretty = TRUE, auto_unbox = TRUE, digits = 8)

cat("Analysis completed\n")
cat(toJSON(summary_out, pretty = TRUE, auto_unbox = TRUE, digits = 5))